#!/usr/bin/perl
# 
# Exercise 9.6
# 
# Write a subroutine that returns a restriction digest, the fragments of DNA left after 
# performing a restriction reaction. Remember to take into account the location of the 
# cut site. (This requires you to parse the REBASE bionet in a different manner. You 
# may, if you wish, ignore restriction enzymes that are not given with a ^ indicating a 
# cut site.)
#
# Answer to Exercise 9.6

# We'll use and modify the answer to Exercise 9.3

use strict;
use warnings;
use BeginPerlBioinfo;     # see Chapter 6 about this module

# Declare and initialize variables
my %rebase_hash = ();

unless(dbmopen(%rebase_hash, 'DBNAME', 0644)) {

	print "Cannot open DBM file DBNAME with mode 0644\n";
}

my @file_data = ();
my $query = '';
my $dna = '';
my $recognition_site = '';
my $regexp = '';
my @locations = ();
my @digest = ();

# If there is a command-line argument, assume it's DNA
if(@ARGV) {
    $dna = $ARGV[0];

# Otherwise, prompt for a FASTA file to open
}else{

    print "Input a FASTA filename: ";
    my $filename = <STDIN>;
    chomp $filename;

    unless(-e $filename) {
    	print "$filename does not exist!\n";
	exit;
    }
    # Read in the file 
    @file_data = get_file_data($filename);

    # Extract the DNA sequence data from the contents of the FASTA file
    $dna = extract_sequence_from_fasta_data(@file_data);
}

# Get the REBASE data into a hash, from file "bionet.110"
%rebase_hash = parseREBASE('bionet.110');

# Prompt user for restriction enzyme names, create restriction map
do {
    print "Search for what restriction site (or quit)?: ";
    
    $query = <STDIN>;

    chomp $query;

    # Exit if empty query
    if ($query =~ /^\s*$/ ) {

        exit;
    }

    # Perform the search in the DNA sequence
    if ( exists $rebase_hash{$query} ) {

        ($recognition_site, $regexp) = split ( " ", $rebase_hash{$query});

	# Find the offset of the cut site in the recognition site
	my $cut_site_offset = index($recognition_site, '^');

        # Create the restriction map
        @locations = match_positions($regexp, $dna);

	# Add the cut site offset to the @locations
	@locations = grep($_ += $cut_site_offset, @locations);

        # Calculate and report the restriction digest to the user
        if (@locations) {
            print "Searching for $query $recognition_site $regexp\n";
            print "A restriction digest for $query cutting at locations:\n";
	    print join(" ", @locations), "\n";

	    # Extract the DNA fragments between the cut sites: the restriction digest
	    for ( my $i = 1, my $j = shift(@locations) ; @locations ; $i = $j, $j = shift(@locations) ) {
		push(@digest, substr($dna, $i-1, $j-$i));
	    }
	    print "The resulting restriction digest:\n";
            print join("\n", @digest), "\n";

        } else {
            print "A restriction enzyme $query is not in the DNA:\n";
        }
    }
    print "\n";
} until ( $query =~ /quit/ );

exit;

